function [dAdv,dzetadalpha,dzetadv,dzetads] ...
    = gradzeta(dPds,dPdv,H,impossible,K,Nu,P,phi,tau,theta0)
% Gradient of stationary preference point distribution (zeta)

alpha = theta0(end-1);


% Calculate transition matrix (same code as stationary distribution calc.)
Pnot0 = cell(1,tau);
for idx1 = 1:tau
    Pnot0{idx1} = P{idx1}(2:end,:,:);
    Pnot0{idx1}(repmat(impossible(2:end,:),[1,1,Nu])) = 0;
end

A = cell(1,tau);
for idx1 = 1:tau
    for idx2 = 1:Nu
        
        a = zeros(K,K);
                
        for idx3 = 1:K

            a(idx3,min(idx3+1,K)) = P{idx1}(1,idx3,idx2) + (1 - phi(2:end-1,idx3)')...
                *Pnot0{idx1}(1:end-1,idx3,idx2) + P{idx1}(end,idx3,idx2)*(idx3==K);
            a(idx3,1) = phi(2:end-1,idx3)'*Pnot0{idx1}(1:end-1,idx3,idx2) ...
                + P{idx1}(end,idx3,idx2)*(idx3==1);            
            if idx3 > 1 && idx3 < K
                a(idx3,idx3) = P{idx1}(end,idx3,idx2); 
            end
            
        end
        
        A{idx1}(:,:,idx2) = a;
        
    end  
end


% Derivative of A wrt baseline utility (v_itk)
dAdv = cell(1,tau);
for idx1 = 1:tau
    dAdv{idx1} = zeros(K,K,H-1,Nu);
    for idx2 = 1:Nu
        for k = 1:H-1
            for idx3 = 1:K
                                
                dAdv{idx1}(idx3,min(idx3+1,K),k,idx2) = dPdv{idx1}(1,k,idx3,idx2)...
                    + (1 - phi(2:end-1,idx3)')*((1-impossible(2:end-1,idx3)).*dPdv{idx1}(2:end-1,k,idx3,idx2))...
                    + dPdv{idx1}(end,k,idx3,idx2)*(idx3==K);
                dAdv{idx1}(idx3,1,k,idx2) = phi(2:end-1,idx3)'*dPdv{idx1}(2:end-1,k,idx3,idx2)...
                    + dPdv{idx1}(end,k,idx3,idx2)*(idx3==1);
                if idx3 > 1 && idx3 < K
                    dAdv{idx1}(idx3,idx3,k,idx2) = dPdv{idx1}(end,k,idx3,idx2);
                end
            end
        end
    end
end


% Derivative of A wrt exponential distribution parameter (s)
dAds = cell(1,tau);
for idx1 = 1:tau
    dAds{idx1} = zeros(K,K,Nu);
    for idx2 = 1:Nu
        for idx3 = 1:K
            
            dAds{idx1}(idx3,min(idx3+1,K),idx2) = dPds{idx1}(1,idx3,idx2)...
                + sum((1 - phi(1:end-1,idx3)).*((1-impossible(1:end-1,idx3)).*dPds{idx1}(1:end-1,idx3,idx2)))...
                + dPds{idx1}(end,idx3,idx2)*(idx3==K);
            dAds{idx1}(idx3,1,idx2) = phi(1:end-1,idx3)'*dPds{idx1}(1:end-1,idx3,idx2)...
                + dPds{idx1}(end,idx3,idx2)*(idx3==1);
            if idx3 > 1 && idx3 < K
                dAds{idx1}(idx3,idx3,idx2) = dPds{idx1}(end,idx3,idx2); 
            end
            
        end
    end
end


% Derivative of zeta wrt baseline utility (v_itk)
dzetadv = cell(1,tau);
b = [alpha; zeros(K-1,1)];
for idx1 = 1:tau
    dzetadv{idx1} = zeros(K,H-1,Nu);
    for idx2 = 1:Nu
        for k = 1:H-1
            dzetadv{idx1}(:,k,idx2) = inv(eye(K) - (1 - alpha)*A{idx1}(:,:,idx2)')...
                *((1 - alpha)*dAdv{idx1}(:,:,k,idx2)')*inv(eye(K) - (1 - alpha)*A{idx1}(:,:,idx2)')*b;
        end
    end
end


% Derivative of zeta wrt exponential distribution parameter (s)
dzetads = cell(1,tau);
for idx1 = 1:tau
    for idx2 = 1:Nu
        dzetads{idx1}(:,:,idx2) = -inv(eye(K) - (1 - alpha)*A{idx1}(:,:,idx2)')...
            *(-(1 - alpha)*dAds{idx1}(:,:,idx2)')*inv(eye(K) - (1 - alpha)*A{idx1}(:,:,idx2)')*b;
    end
end


% Derivative of zeta wrt attrition rate (alpha)
dzetadalpha = cell(1,tau);
dbdalpha = [1; zeros(K-1,1)];
for idx1 = 1:tau
    for idx2 = 1:Nu
        dzetadalpha{idx1}(:,:,idx2) = -inv(eye(K) - (1 - alpha)*A{idx1}(:,:,idx2)')...
            *A{idx1}(:,:,idx2)'*inv(eye(K) - (1 - alpha)*A{idx1}(:,:,idx2)')*b ...
            + inv(eye(K) - (1 - alpha)*A{idx1}(:,:,idx2)')*dbdalpha;
    end
end